These slides and all other course materials can be found at
github.com/brry/rhydro #slides
To get all the material including the datasets and presentation source code, we recommend to download the whole github course repository.
This is an R Markdown Notebook.
For discussions, please visit the Hydrology in R Facebook group.
Before running the code blocks below, we suggest to get package installation instructions by running:
source("https://raw.githubusercontent.com/brry/rhydro/master/checkpc.R")
Aim and contents of this workshop
We want to:
We can not:
We have prepared:
rmarkdown, R notebook)sf + leaflet + mapview + OSMscale)animation + extremeStat)airGR
Before we get started, please let us know your current R knowledge level by filling out the short survey at
bit.ly/knowR
Using R as GIS (reading a rainfall shapefile + Kriging, sf + leaflet + mapview + OSMscale)
Berry Boessenkool
Reading shapefiles with maptools::readShapeSpatial and rgdal::readOGR is obsolete.
Instead, use sf::st_read. sf is on CRAN since oct 2016.
Main reaction when using sf: “Wow, that is fast!”
Download the shapefile or better: download the whole github course repository
rain <- sf::st_read("data/PrecBrandenburg/niederschlag.shp")
Reading layer `niederschlag' from data source `/home/berry/Dropbox/R/rhydro/presentations/data/PrecBrandenburg/niederschlag.shp' using driver `ESRI Shapefile'
converted into: POLYGON
Simple feature collection with 277 features and 1 field
geometry type: POLYGON
dimension: XY
bbox: xmin: 3250175 ymin: 5690642 xmax: 3483631 ymax: 5932731
epsg (SRID): NA
proj4string: +proj=tmerc +lat_0=0 +lon_0=15 +k=0.9996 +x_0=3500000 +y_0=0 +ellps=GRS80 +units=m +no_defs
centroids <- sf::st_centroid(rain)
centroids <- sf::st_coordinates(centroids)
Static plot:
plot(rain[,1])
Static map:
prj <- sf::st_crs(rain)$proj4string
centroids <- as.data.frame(centroids)
cent_ll <- OSMscale::projectPoints(Y,X, data=centroids, to=OSMscale::pll(), from=prj)
map_static <- OSMscale::pointsMap(y,x, cent_ll, fx=0.08, type="maptoolkit-topo", zoom=6)
Downloading map with extend 11.029332, 14.981464, 51.10418, 53.765423 ...
Done. Now plotting...
Interactive map:
library(leaflet)
cent_ll$info <- paste0(sample(letters,nrow(cent_ll),TRUE), ", ", round(cent_ll$x,2),
", ", round(cent_ll$y,2))
leaflet(cent_ll) %>% addTiles() %>% addCircleMarkers(lng=~x, lat=~y, popup=~info)
Interactive map of shapefile:
# devtools::install_github("environmentalinformatics-marburg/mapview", ref = "develop")
library(berryFunctions) # classify, seqPal
col <- seqPal(n=100, colors=c("red","yellow","blue"))[classify(rain$P1)$index]
mapview::mapview(rain, col.regions=col)
Plot original points colored by third dimension:
pcol <- colorRampPalette(c("red","yellow","blue"))(50)
x <- centroids$X # use cent_ll$x for projected data
y <- centroids$Y
berryFunctions::colPoints(x, y, rain$P1, add=FALSE, col=pcol)
Calculate the variogram and fit a semivariance curve
library(geoR)
--------------------------------------------------------------
Analysis of Geostatistical Data
For an Introduction to geoR go to http://www.leg.ufpr.br/geoR
geoR version 1.7-5.2 (built on 2016-05-02) is now loaded
--------------------------------------------------------------
geoprec <- as.geodata(cbind(x,y,rain$P1))
vario <- variog(geoprec, max.dist=130000) # other maxdist for lat-lon data
variog: computing omnidirectional variogram
fit <- variofit(vario)
variofit: covariance model used is matern
variofit: weights used: npairs
variofit: minimisation function used: optim
initial values not provided - running the default search
variofit: searching for best initial value ... selected values:
sigmasq phi tausq kappa
initial.value "1325.81" "19999.05" "0" "0.5"
status "est" "est" "est" "fix"
loss value: 104819060.429841
plot(vario)
lines(fit)
Determine a useful resolution (keep in mind that computing time rises exponentially with grid size)
# distance to closest other point:
d <- sapply(1:length(x), function(i)
min(berryFunctions::distance(x[i], y[i], x[-i], y[-i])) )
# for lat-long data use (2017-Apr only available in github version of OSMscale)
# d <- OSMscale::maxEarthDist(y,x, data=cent_ll, fun=min)
hist(d/1000, breaks=20, main="distance to closest gauge [km]")
mean(d/1000) # 8 km
[1] 8.165713
Perform kriging on a grid with that resolution
res <- 1000 # 1 km, since stations are 8 km apart on average
grid <- expand.grid(seq(min(x),max(x),res),
seq(min(y),max(y),res))
krico <- krige.control(type.krige="OK", obj.model=fit)
#krobj <- krige.conv(geoprec, loc=grid, krige=krico)
#save(krobj, file="data/krobj.Rdata")
load("data/krobj.Rdata") # line above is too slow for recreation each time
Plot the interpolated values with or an equivalent (see Rclick 4.15) and add contour lines.
par(mar=c(0,3,0,3))
geoR:::image.kriging(krobj, col=pcol)
colPoints(x, y, rain$P1, col=pcol, legargs=list(horiz=F, title="Prec",y1=0.1,x1=0.9))
points(x,y)
plot(rain, col=NA, add=TRUE)
River discharge time-series visualisation and extreme value statistics (animation + extremeStat)
Berry Boessenkool
Exploratory Data Analysis including flow duration curve and trend analysis on time-series
Shaun Harrigan
Please give us feedback at bit.ly/feedbackR
For discussions, please use the Hydrology in R Facebook group.